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Abstract 

Implicit-explicit (IMEX) time stepping methods can efficiently solve differential 
equations with both stiff and nonstiff components. IMEX Runge-Kutta methods and 
IMEX linear multistep methods have been studied in the literature. In this paper we 
study new implicit-explicit methods of general linear type (IMEX-GLMs). Wc develop 
an order conditions theory for high stage order partitioned GLMs that share the same 
abscissae, and show that no additional coupling order conditions are needed. Conse- 
quently, GLMs offer an excellent framework for the construction of multi-method inte- 
gration algorithms. Next, we propose a family of IMEX schemes based on diagonally- 
implicit multi-stage integration methods and construct practical schemes of order three. 
Numerical results confirm the theoretical findings. 
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1 Introduction 



Implicit-explicit (IMEX) time integration schemes are becoming increasingly popular for 
solving multipliysics problems with both stiff and nonstiff components, which arise in many 
application areas such as mechanical and chemical engineering, astrophysics, meteorology 
and oceanography, and environmental science. Examples of multiphysics problems with both 
stiff and nonstiff components include advection-diffusion-reaction equations, fluid-structure 
interactions, and Navier-Stokes equations. Such problems can be expressed concisely as the 
system of ordinary differential equations (ODEs) 

y' = fit^y) + git^y) to<t<tF, y(to) = i/o, (i) 

where / corresponds to the nonstiff term, and g corresponds to the stiff term. In case of sys- 
tems of partial differential equations (PDEs) the system ([T]) appears after semi-discretization 
in space. 

An IMEX scheme treats the nonstiff term explicitly and the stiff term implicitly, therefore 
combining the low cost of explicit methods with the favorable stability properties of implicit 
methods. IMEX linear multistep methods have been developed in [H [21 [3], and IMEX 
Runge-Kutta methods have been built in [H (Sj [6l [7] . 

The general linear method (GLM) family proposed by J.C Butcher [S] generalizes both 
Runge-Kutta and linear multistep methods. The added complexity improves the flexibility 
to develop methods with better stability and accuracy properties. While Runge-Kutta and 
linear multistep methods are special cases of GLMs, the framework allows for the construc- 
tion of many other methods as well. Here we focus on the diagonally implicit multistage 
integration methods (DIMSIM) [9J, which are both efficient and accurate, and great poten- 
tials for practical use. GLM can overcome the limitations of both linear multistep methods 
(lack of A-stability at high orders) and of Runge-Kutta methods (low stage order f leading 
to order reduction). A complete treatment of GLMs can be found in the book of Jackiewicz 
[lOj. 

In this study we develop the concept of partitioned DIMSIM methods, and develop an 
order conditions theory for a family of such methods. This shows that partitioned GLM is a 
great framework for developing multi- methods. Next, we propose a new family of implicit- 
explicit methods based on pairs of DIMSIMs, and develop second and third order methods 
on this class. 

In our earlier work [TT] we have developed second order IMEX-GLM schemes. While this 
paper was under study, we became aware of an effort to construct IMEX-GLM schemes for 
Hamiltonian systems |12j . 

The paper is organized as follows. Section |2] reviews the class of general linear methods. 
The new concept of partitioned DIMSIM schemes is proposed in Section [3| and the order 
conditions theory is developed. IMEX-DIMSIM schemes are constructed in Section |4j Linear 



stability is analizes in section 4.5, and Prothero- Robinson convergence in section [475| IMEX 



methods of second and third order are built in Sections |5.1| and |5.2[ respectively. Numerical 
results for van der Pol system and for the two dimensional gravity waves equations are 
presented in Section [6j Section [7] draws conclusions and points to future work. 
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2 General linear methods 



2.1 Representation of general linear methods 

Consider the initial value problem for an autonomous system of differential equations in the 
form 

y'{t) = f{y), to<t<tF, y(to)=2/o, (2) 

with f : ^ W'' and y{t) E R'^. GLMs [TO] for ^ can be represented by the abscissa 
vector c e M*, and four coefficient matrices A G W^", U G W^"^, B G W^" and V G 
which can be represented compactly in the following tableau 



rxr 



A 


U 


B 


V 



On the uniform grid t„ = to + ^h, n = 0,1, . . . , N, Nh = tp — to, one step of the GLM reads 

s r 

Yi = h^aijf(Yj) + ^Ui^jyf'^\i = l,...,s, (3a) 

s r 

= hY,h,f{YJ) + J2v^,yf''\^ = ^,■■■,r, (3b) 
i=i i=i 

where s is the number of internal stages and r is the number of external stages. Here, h is 
the step size, Yi is an approximation to y(t„-i + Cih) and |/|"' is an approximation to the 
linear combination of the derivatives of y at the point tn- The method ^ can be represented 
in vector form 

Y = h{A® Irfxd) F{Y) + (U ® Irfxd) (4a) 
= (B®I,xd) F(r) + (V®I,xd) 2/'"-'], (4b) 

where I^xd is an identity matrix of the dimension of the ODE system. 
2.2 Stability considerations 

The linear stability of method ^ is analyzed in terms of its stability matrix 

M{z) =V + zB (I,x. - zA)-' U , (5) 
and the corresponding stability function 

p{w, z) = det{wlrxr — ^{z)), (6) 



where w,z EC A desirable property is the inherited Runge-Kutta stability [T3|,[T4]. This 
means that the stability function ^ has the form 

p{w, z) = w'-^ {w - R{z)) , (7) 

where R{z) is the stability function of Runge Kutta method of order p = s. 
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2.3 Accuracy considerations 

We assume that the components of the input vector yY ^' for the next step in ([3]) satisfy 



Y,%kh''y^''\tn-i) + 0{h^+^), z = 1, . . . ,r, 



fc=0 



for some real parameters (jj^fc, z = 1, . . . , r, = 0, 1, . . . ,p. 

The method ^ has order p if the output vector ?/|"' satisfies 



= Y.'i^,khV'\tn) + ^ = 1, . . . ,r, 

fc=0 



(9) 



for the same parameters qi^k of ([8j). 

The method ^ has stog'e order q if the internal stage vectors y/"' are approximations of 
order q to the solution at the time points tn-i + Cih 



= y(t„„i + cih) + o{h''+'), i = i,...,s. 

We collect the parameters qi k in the matrix W for convenience 



(10) 



W = [qo qi 



qp] 





91,1 ■ ■ 




Q2fi 


12,1 ■ ■ 


■ l2,p 









Theorem 1 (GLM order conditions [lO]). Assume that y'" -"^^ satisfies Then the GLM 
has order p ^ and stage order q = p [T^ if and only if 



e'w{z) = zBe^' + Yw{z) + 0{z^+^) 



(12a) 
(12b) 



where 



j=0 



For stage order q = p — 1 condition ( |12a[ ) is replaced by 

cP 



zKe""- + \}w(z) + 



Ac 



jo! [p — \)\ 

Proof. See Butcher [S] and Jackiewicz [IQi Section 2.4]. 



Uqp + 0(;2P+i) . 



(12c) 

□ 



It is shown in [10] that a GLM ([S]) has order j9 and stage order q with g = p = r = s if 
and only if 

B = Bo - ABi - VB2 + VA, (13) 
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where the matrices Bq, Bi, B2 € M*^* are defined by 



with 

s 

2.4 Starting and ending procedures 

Assumption (|8| requires to compute the initial vector yt*^' by a starting procedure satisfying 

2/r'=E9.fc/^¥'Hio) + 0(/^^+'), ^ = l,...,r. (15) 

fc=0 

Dense output is based on derivative approximations of the form 

s r 

h'y^'\tn) ^ + k = 0,l,...,r. (16) 

i=0 i=0 



It is shown in [H [15] that (16) is accurate within (9(/i^+^) if and only if 



[1, 2, ... , zPfe' = zBe^' + Vw{z) + 0{zP+^) (17) 



where B = j] and V = [7^,,]. The termination procedure uses (16) with A; = to generate 
the solution at the last step 

s r 

y{tn) ^ J^MXYi) + T.^oJr'^. (18) 

i=0 j=0 

2.5 Diagonally implicit multistage integration methods 

Diagonally implicit multistage integration methods (DIMSIMs) are a subclass of GLMs char- 
acterized by the following properties [8]: 

1. A is lower triangular with the same element Oj^j = A on the diagonal; 

2. V is a rank-1 matrix with the nonzero eigenvalue equal to one to guarantee preconsis- 
tency; 

3. The order p, stage order q, number of external stages r, and number of internal stages 
s are related by g G {p — l,p} and r E {s,s + 1}. 

In this work we focus on DIMSIMs with p = q = r = s,\J = Isxs, and V = Isv'^, where 
v'^ Is = I [lOj. DIMSIMs can be categorized into four types according to [8]. Type 1 or type 
2 methods have aij = for j > i and are suitable for a sequential computing environment, 
while type 2 and type 3 methods have j = for j ^ i and are suitable for parallel 
computation. Methods of type 1 and 3 are explicit (oj^j = 0), while methods of type 2 and 
4 are implicit (aj^j = A 7^ 0) and potentially useful for stiff systems. 
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3 Partitioned general linear methods 

Consider the partitioned system of ODEs 

f{i}{y{i}, - ■ ■ ,y{N}) 
f{N}{y{i},---,y{N}) 



(19) 



where the solution vector is separated into components y{m}, ^ = ■ ■ ■ ,N, each of which 
may be itself a vector. 

A partitioned general linear method solves (19) by applying a different GLM to each 
component. If not explicitly stated otherwise , we use the subscript {m} to denote the 
coefficients specific to the m-th component of the partitioned system. We have the following 

Definition 1 (Partitioned GLM). One step of a partitioned GLM has the form 



Y, 



{m}i 



h ^ a{m}ij f{m}{y{l}j,y{2}j, • • • , y{N}j) 



(20a) 



\ ^ [n-l] 



Z = 1,...,S, 



y{m}i 



h ^ b{m}ij f{m}(y{i}j, y{2}j, • • • , y{N}j] 

r 



(20b) 



where a^m}i,j,U{m}ij,b{m}i,j, o-nd c^m}i for m = 1, . . . , N represent the coefficients of N dif- 
ferent GLMs. 



Definition 2 (Internal consistency). A partitioned GLM (20) is internally consistent if all 
component methods share the same abscissae, C|m}j = Cj for m = 1, . . . , N . 

Internal consistency means that all stage components approximate the solution compo- 
nents at the same time point, i.e., [^{iij, , • • • , ^{Af}i] ~ y(tn + Cjh), for all j = 1, . . . , s. An 
internally consistent partitioned GLM method (20) can be represented compactly as 



-{m} 





U{m} 


B{.m} 





(21) 



Definition 3 (Order of partitioned GLM). Assume that each component of the input vector 
satisfies (Isj) 



<YJ = E ^{-K'^^' yil}('n-i) + Oih^^'), z = 1, . . . , r . 



(22) 



fc=0 
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The partitioned GLM (20) has order p if each component of the output vector satisfies 



m 



(23) 



k=0 



for the same parameters q{m}i,k o,s in (22). The partitioned GLM (20) has stage order q if 
each component of the internal stages 1^'"' satisfies 



'{^}^ = y{n.}{tn-l + C{^}^ h) + 0{h''+^), I 



1 S . 



m 



1 N. 



(24) 

Theorem 2 (Order conditions for partitioned GLMs). Assume that each component y^lln}j 
satisfies Q). Then the internally consistent partitioned GLM (21) has order p (23) and stage 
order q ^ {p — l,p} {24) if and only if each component method (Aj^}, Bj^}, U{m}, V{m}) 
has order p ^ and stage order q (10). 

Remark 1. Each component method needs to independently meet its own order conditions 
(12). No additional "coupling" conditions are needed for the partitioned GLM (i.e., no order 
conditions contain coefficients from multiple component schemes). 



Proof. We first prove the "only if" part: if the partitioned GLM satisfies (23)-(24) with order 
p stage order q & {p — l,p}, then each component method satisfies its own order conditions 
with the same p and q. This can be seen immediately by employing the same 
component method for all partitions, {Ai^k},'B{k},V{k},^{k}) = (Aj^}, B{^}, U{m}, Vj^}) 
for k = 1, . . . , N. The partitioned method (21 ) is the traditional GLM method (Aj^} , B{m}, 
U{m}7 ^{m}) and has to satisfy the traditional order conditions ^ and (10). 

We next prove the "if" part: if each component method satisfies ([9])-(10) with order p 
stage order q E {p — l,p}, then the partitioned GLM (21) has order p and stage order q. 
Denote 



and 



y{tn-i + Cjh) 



Y{N}j_ 

?/{i}(t„_i + Cjh) 

y{N}{tn-l + Cjh) 



Y 



y{tn-i + ch) 



y{tn-i + cih) 



y{tn~i + Csh) 



Consider the stage equations of the individual method m with exact solution arguments 

s 

y(tn-i + Cih) = /i ^ a{rn}i,j f {y(tn-i + Cjh)) (25) 

' p \ 

Yl Q{m}^,kh''y^'\tr,.l) + O {h'^+') , t = 1, . . . , S . 



i=i 




fc=0 



The error size is given by the stage order q of each individual method (10). Using the 
assumption (22) each component of the sum ^^=0 ^{™}*.^^'^^{m}(^"-i) "^^^ replaced by 
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the numerical approximations which differ from their exact counterparts by 0{h^~^^)] 

therefore their use in (25) does not change the asymptotical error size. The m-th component 



of relation (25) then reads 



y{m}{tn-l + Cih) = Ct{m}i,j f{m} {y{tn-l + Cjh)) 



(26) 



r 



Subtracting (26) from the stage equation (20a) gives 



y{m)i - y{m}{tn~l + Cih) = a{m}ij {f{m}iYj) - /{m}(?/(t„-l + Cjh))) + 0{h''^^) 

and therefore 

\\Y{m}i - y{m}{tn-l + Q^)|L - ^ ll^wlL ^™ ~ y{tn-l + ch)\\^ + C>(/l'?+^) 

where Lm is the Lipschitz constant of f{m}- It follows that [lO] 

\\Y - y{tn., + ch)\\^ = 0{h'^+') (27) 
for all sufficiently small step sizes 



h < T = { max 1 1 A 



V m 



"{mil loo 



Equation (27) proves the stage order of the partitioned GLM method. 



Continuing, (27) implies that 



hf{m}{Y,) = hf{m}{y{tn-l+C,h)) + 0{h'l+^). 

= hf{m}{y{tn-l+Cih)) + 0{K'+^) 



(28) 



where we have used the fact that g + 2 >p + l. Consider the solution step of the individual 
method m with exact solution arguments 



fe=0 



^q{m}i,kh^y^''\tn) = hJ2hm}i,jfiyiin-i + Cjh)) (29) 

r / p \ 



,k=0 



for 2 = 1, . . . , r, where the size of the error term reflects the fact that each individual method 



has order p. 
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Use (28) and the assumption (22) into the m-th component of (29) to obtain 

p s 



(30) 



fc=0 



y';iy^ + 0{h^'-') ^ = l,...,r, 



(31) 



The last equahty follows from the partitioned method solution equation (20b). This estab- 
lishes the order p of the partitioned GLM. 

□ 



4 Implicit-explicit general linear methods 
4.1 Construction procedure 

The derivation of IMEX-GLM schemes relies on the partitioned GLM theory developed in 
Section [3| W e transform the additively partitioned system ([T]) into a component partitioned 
system ( 19 ) via the following transformation jl] 



y = x^+z, 

x' = f{x,z) = f{x + z), 
z' = g{x, z) = g{x + z) . 



Equation (32a) is discretized with an explicit (type 1) GLM 



i-l 



Xi = aij f{Xj + Zj) + Uij xf ^\ i = l, 



i=i 

s 



r 



X- 



[n-l\ 



I = 1, 



Similarly, equation (32b) is discretized with an implicit (type 2) GLM 



(32a) 
(32b) 



(33a) 
(33b) 



[n-l] 



Zi = h^aij g{Xj + Zj) + ^ Uij z^ "\ i = 1, . . . , s 

zf = hJ2h 9{X, + Z,) + 4""" 



i = l,...,r. 



(34a) 
(34b) 
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Combining (33) and (34) we obtain 

'i-l 



Xi + Zi 



(35a) 



in-l] 



[n-1] 



Z = 1, . . . ,S, 



h J2 fix, + z,) + 9{X, + z, 



(35b) 



r 



[n.-l] , ^ [n-1] 



i = l,...,r, 



We consider pairs of explicit ( 33 ) and implicit ( 34 ) schemes that 



• share the same abscissa vector c = c such that the partitioned GLM is internally 
consistent, and 

• share the same coefficient matrices U = U and V = V. 

For this class of schemes all internal stage vectors can be combined. Specifically, let Yi = 



Xi + Zi and yi = Xi + Zi. The scheme (35) becomes the following method/ 



Definition 4 (IMEX-GLM methods). One step of an implicit- explicit general linear method 
applied to ([T]) advances the solution using 



Yi = h^aijf{Yj) + hJ2^i,j9{Yj) + ^Uijy^- 



[n-1] 



i = l,...,s, 



i=i 

s 



[n-1] 



i = l,...,r. 



(36a) 
(36b) 



We note that in (36) and need not to be known individually once they are 
initialized ine the ffist step. The combined solution ?/|"^ = x^"' + is advanced at each step 



as regular GLMs do. The IMEX-GLM (36) is represented compactly by the Butcher tableau 

c A A U 

b | b | v ■ (3^) 

4.2 Starting procedures 



An IMEX GLM (|36|) of order p requires a starting procedure that approximates linear com- 

r 

f^ = J2^^,kh'z'^'\to) + 0{hn (38) 



binations of derivatives as follows 

r 

J2%kh''x^''\to) + 0{hP) and 
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fe=0 



fc=0 



respectively, where 



i-1 



Ac 



i! 



Ac 



(39) 



Thus 



x{to) + z{to) + qi,ihx'{to) + %ihz'{to) 

r r 

+ J2 qi,kh'x^''\to) + J] qi,kh'z^''\to) 

k=2 k=2 

yo + qi,ihf{yo) + qi,ihg{yo) 

r r 



k=2 



k=2 



Evaluation of the first three terms is straightforward. But approximations of the other 
terms containing derivatives x^''\to) and y^''\to) for A; > 2 requires additional work if their 
analytical expressions are difficult to obtain. 

To initialize an IMEX GLM we approximate independently the vectors h''x^''\to), h^z^^\to), 
k — 1, . . . , r, using finite differences and the solution information provided by several steps 
of an IMEX Runge-Kutta method. 

For better accuracy, the IMEX RK method uses a small step size r < h, and produces 
the numerical solutions y^*^"^* ^ y{tQ + ir). In the following we show how to compute the 
terms t'^ x'^'^^to)] each of these terms is then rescaled by {h/rY to reflect the integration step 
h. We have that 



Tx'ito) 
T^x"{to) 




■ x'ito) ' 
x'ih) 


+ 0{t'+^) = tD 


f{yo) 
fiyf"'') 


+ 0{t'+') (40) 


T'-xW(to) 




X'itr) 









where the coefficient matrix D G W^"^ is derived by expanding the right hand side in Taylor 
series and comparing the coefficients of each term. For the cases r = 2 and r — 3 the 
coefficients are 



D 



(r=2) 



1 

-1 1 



and D 



(r=3) 



1 

-3/2 
1 




2 





■1/2 
1 



respectively. The same procedure is applied to obtain r'^z^'^^to). We note that the initializa- 
tion procedure requires the function values f{y) and g{y) evaluated at the starting solution 
steps yf^^'^^^ and that there is no need to compute Xi or Zi separately. 
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4.3 Termination procedures 



To generate the solution at the last time step yitp) using ( 18 ) a general termination procedure 
reads 

s r 

y{tn) ^ Y.Miy^ + Y.^o,xf-'^ (41a) 

i=0 j=0 



1=0 j=0 



In order to avoid separate evaluations of x^" ^' and ^j" we require that 7oj- = 7o,j for all 
j. In this case the termination procedure reads 

s s r 

y{tn) ^ J2Miy^ + J2My^ + J2^oWr'^. (4ib) 

i=0 1=0 j=0 

For explicit (type 1) GLMs, choosing the first abscissa coordinate ci = implies that 
qifi = 1 and qi^i = for z > 1 due to order conditions. The first element of the output vector 
is exactly the solution at the current step, i/^"' ~ y{tn)- In this case, /3o is equal to the first 
row of the coefficient matrix B, and 70 is the first row of V. 

For implicit (type 2) GLMs, there are usually sufficiently many free parameters in B and 



V that remain after satisfying (13). These free parameters could be chosen in such a way 



that the implicit GLM shares the same coefficients 70 with the explicit GLM. The difficulty 
of computing terms a^^""^^ and -z]""^' individually can therefore be avoided. 

4.4 Linear stability analysis 



For convenience, we write the IMEX-GLM (36) in the vector form 



Y = hAF{Y) + hAG{Y) + Vy^''-^^ (42a) 

yl"] = hBF{Y) + hBG{Y)+yy^''-^K (42b) 

We consider the generalized linear test equation 

y' = ^y + h, t>o, (43) 

where ^ and ^ are complex numbers. We consider C,y to be the nonstiff term and the stiff 
term, and denote w = hS, and w = hC,. 

Applying (42) to the test equation (43) leads to 

Y = /i (^^A + ^a) F + Uy["-^], (44a) 

= /i f^B + fs) r + Vyl""^] . (44b) 



Assuming I^xs ~ wA — wA is nonsingular we obtain 
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where the stabihty matrix is defined by 

M{w,w) = V + (^wB + wB^ (isxs-wA-wA^ ^U. (45) 

Let 5 C C and 5 C C be the stabihty regions of the exphcit GLM and of the imphcit GLM, 
respectively. The combined stability region is defined by 

l^weS,weS : p{M{w, w)) C S xS cCxC. (46) 

For a practical analysis of stability we define a desired stiff stability region, e.g., 

Sa = {weSnC' : |Im(w)| < tan(a) |Re(w)|}, 
and compute numerically the corresponding non-stiff stability region: 

Sa = e S : p{M{w,w)) <1, VwG5a}. (47) 

The IMEX-GLM method is stable if the constrained non-stiff stability region Sa is non-trivial 
(has a non-empty interior) and is sufficiently large for a prescribed (problem-dependent) value 
of a, e.g., a = 90°. 



4.5 Prothero-Robinson convergence 

We now study the possible order reduction for very stiff systems. We consider the Prothero- 
Robinson (PR) [16j test problem written as a split system ([T| 

y' = fi{y-<P{t))+<P'it) , /i<0, y(O) = 0(O), (48) 

9{y) f(y) 

where the exact solution is y{t) = (j){t). A numerical method is said to be PR-convergent 



with order p if its application to (48) gives a solution whose the global error decreases as 
0{hP) for h and hfi -oo. 

Theorem 3 (Prothero-Robinson convergence of IMEX-GLM). Consider the IMEX GLM 



method (36). Without loss of generality we consider that U = I. The explicit part is of 



order p and stage order q & {p — and the implicit part has order p = p and stage order 



g G {j5 — Assume that hfi E S for all h > 0. Then the IMEX GLM method (36) is 

PR-convergent with order min(p, q). 

Remark 2. // the explicit stage order is q = p, then the PR order of convergence is p. It 



is convenient to construct IMEX GLM methods (36) with explicit stage order q = p, even if 



q = p — 1, as such methods do not suffer from stiff order reduction on the PR problem. 
Proof. Let 

= (t„_i + ch) = [0(t„_i + Ci /l), . . . , 0(t„_i + c, h)f . 

and 

= [(/)(t„_i),/i0'(t„_i),...,/i^'0(^')(t„_i)]^ . 
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The method (36) apphed to (48) reads: 

r M = hA + hiiA (r ["1 - + U yl"^^] , (49a) 
= /iB0'["l + /i/iB (yM + Vy["-i] . (49b) 

Consider the global stage errors 

To obtain the global error in we consider separately the global errors in the nonstiff and 
stiff components: 

gnonstiflf _ ^[n] _ ^ h'' X^''^ (t„) , 

= 0N_^H_ (#•) - x^'^)) (t„) 



since the exact solution of the nonstiff system is x{t) = (j){t). Consequently, the total error is 



3 nonstiff _i_ stiff 



Write the stage equation (49a) in terms of the exact solution and global errors 

p 

fc=0 

to obtain 

(i,xs-/^/^a) = e„_i + /iA0'(t„_i + c/i) (50) 

p 

+U ^ qfc <P^'\tn-l) - <p{tn-l + Ch) . 
k=0 

The exact solution is expanded in Taylor series about t„_i: 

(f){tn-l + ch)-ls(f){tn-l) = I] -^^^''^in-l) , 

fc=l 

/.0'(t._i + c/i) = ^^^^0W(t„_,). 

k=l 
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Inserting the above Taylor expansions in (50) leads to 



+ J2{kA c'^-i + fc! U qfc - c^) — {t 



k=l 



n—1 ) 



where q is the stage order of the explicit method. We have used the facts that qo = l^, 
U Is = Is, and the order conditions (12a) and (12c) for the cases where g = p and q = p—1, 
respectively. 

Similarly, we write the solution equation (49b) in terms of the exact solution and global 
errors: 

p 

er, + J2^kh'<P^''\Q = /iB0'(t„_i + c/i) + /i/iBEW+Ve[„_i] 

fc=0 



+V^qfc/l'^0W(t„_ 
After rearranging the expression we obtain 

p 

+h B + c/i) + V ^ qfc [tn-i] 

k=0 



k=0 



By Taylor series expansion we have 



k=0 



k=o \e=o 



and therefore 



M(V) e 



n-l 



A;Bc'^-^ + A;!Vqfc-A;!^ 



A:=l 



Qk-e \ 
i\ / k\ 



(tn-l) + O (/l^+l) 



The order condition ( 12b ) of the nonstiff scheme reads 



e'w{z) = zB e'' + Yw{z) + O {zP+^) 



e>o k=o 



k=0 



k\ 



k=0 



(51) 
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Identification of powers of leads to 



B 



{k-l)\ 



+ VcikZ^ , k = l,. 



Tfie error recurrence (51) becomes 



(52) 



Assume tliat the initial error is Cq = 0{h^). The error amplification matrix lSA{hfi) is the 
stability matrix of the implicit method. Therefore its spectral radius is uniformly bounded 
below one for all argument values hn of interest. By standard numerical ODE arguments [IT] 
the equation (52) implies convergence of global errors to zero at a rate ||e„|| = O (^/i^^^P'?)). 

□ 



Construction of implicit-explicit methods of orders 
two and three 



We now construct IMEX-DIMSIM methods as summarized in Section 2.5 

r = s, U = Isxs, and V = 1., w where f ^ 



focus on DIMSIMs with p = q 



Specifically, we 

1. = 1 m. 



5.1 Two-stage, second-order pairs with p = q = r = s = 2 

The pair of explicit and implicit schemes developed in [llj is named IMEX-DIMSIM-2A and 
consists of a type 2 DIMSIM from [8J with the same stability of SDIRK method of order 
2, and a type 1 derived DIMSIM. Both of them share the same abscissa vector c = [0, 1]^ 
and the same coefficient matrix V. The IMEX-DIMSIM-2A coefficients in the tableau (37) 
representation are 












2-V2 
2 





1 





1 


2 





2V2+6 
7 


2-v^ 
2 





1 




3v^-l 
4 

3\/2-3 
4 


3-V2 
4 

4 


73-34^2 
28 

87-48^2 
28 


4^2-5 
4 

-45+34^2 
28 


3^2-3 
4 

3-v^ 
2 


1-V2 

4 

V2-I 
2 



The choice of A = (2 - \/2) /2 ensures the type implicit part of IMEX-DIMSIM-2A is 
L-stable. Inherited Runge-Kutta stability is a desirable property, but there are not enough 
free parameters to enforce this property on both methods of the IMEX pair at the same 
time. 

For a given implicit scheme we construct the explicit method by maximizing the con- 



strained stability region (47). We have observed that simply maximizing the explicit stability 
region S is insufficient and can lead to a very poor constrained stability region for the IMEX 
method. The matrix B can be determined by A, c and V according to the order condition 



( 13 ). The only free parameter is 02,1 in matrix A, and it is chosen such as to maximize IMEX 
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2p 
1.5 
1 - 
0.5 
0- 
-0.5 
-1 - 
-1.5 




(a) Stability region S of the im- (b) Stability region S of the ex- (c) Constrained stability region 
plicit method plicit method Sa (47) for a = 90° 



Figure 1: Stability regions for the IMEX-DIMSIM-2B pair 

stability. First, we use a Matlab Differential Evolution package as a heuristic for global 
optimization to generate a starting point. Then we run the Matlab routine fminsearch 
multiple times until the result converges; each run is initialized with the previous result. 
The resulting stability regions are reported in Figure [T] 

This procedure led to another explicit scheme that maximizes the IMEX stability 




1.5 



B 



V2 3-V2 
2 4 

V2-I 3-V2 
2 4 



U and V are the same. We call the new pair IMEX-DIMSIM-2B. The termination procedure 

bl,2, 70,1 = ^^1,1, 70,2 = Vl,2- 



(41) has the following parameters 

^0,1 



^1,1) /3o,2 

Solving the condition ( [l7| ) gives 

73 - 34^2 43 - 31v^ 



7o,i 



28 ^ 28 
3-V2 2-V2 



9, A 



0,2 



-I + 2V2 -4 + 3v^ 

+ : 9, 



4 



-9, 70,2 



V2-1 V2-2 



The choice of the free parameter g = leads to 70,1 = 70,1, 70,2 = 70,2, and (41b). 



5.2 Three-stage, third-order pairs with p = q = r = s = 3 

We construct two implicit-explicit pairs named IMEX-DIMSIM-3A and IMEX-DIMSIM-3B 
starting from two existing implicit methods. All coefficients are obtained from the numerical 
solution of order conditions using Mathematica. The calculations are performed with 24 
digits of accuracy such as to reduce the impact of roundoff errors on the resulting coefficient 
values. 



^http: / /www. mathworks.com/matlabcentral/fileexchange/18593-differential-evolution 
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2p 
1.5 
1 - 
0.5 
0- 
-0.5 
-1 - 
-1.5 





(a) Stability region S of the im- (b) Stability region S of the ex- (c) Constrained stability region 
plicit method 



plicit method 



Sa (47) for a = 90° 



Figure 2: Stability regions for the IMEX-DIMSIM-3A pair of schemes 



IMEX-DIMSIM-3A. According to [IE] there are five A-stable type 2 DIMSIMs with the 
choice A = 1/2 and c = [0, 1/2, 1]^. We select the implicit component in Table [l] which has 
a balanced set of coefficients. 

The explicit component is obtained by a numerical maximization of the constrained 
stability region, as discussed in the previous section. The resulting coefficients are shown in 
Table |3| The IMEX stability regions are drawn in Figure |2j 

The termination procedure (41) is given by 



hi 
7o,i 



1.01640094894605, /3o,2 = 0.632229903531054, 

^1,1) Po,2 = ^1,2; PqZ = bl3, 

70,1 = ^^1,1, 70,2 = 70,2 =^^1,2 , 70,3 = 703 = ^1,3- 



03 



0.0919425241172364, 



IMEX-DIMSIM-3B. The choice of A = 0.435866521508459 and c = [0, 1/2, 1]^ leads to 
the L-stable type 2 DIMSIM reported in |T8j. The coefficients of the implicit component are 
presented in Table |2j 

The type 1 component is shown in Table |4} The IMEX stability regions are drawn in 
Figure [3j 



The coefficients (3 and 7 of the termination procedure (41) are equal to the first rows of 
matrices B and V, respectively. In addition 



0,1 



0.833790728250125, /3o2 = 0.645998912146314, /3o.3 = 0.120039435995489. 



6 Numerical results 

We test the IMEX-GLM methods on two test problems. The first one is the van der Pol equa- 
tion, a commonly used small ODE system that emphasizes convergence under stiffness. The 
second test is a PDE problem arising in atmospheric modeling. We implemented our algo- 
rithms in a discontinuous Galerkin finite element model developed by Blaise et al. [19], which 
has efficient parallel scalability. We report the results obtained with IMEX-DIMSIM-2B and 
IMEX DIMSIM-3B methods, since they have the better accuracy and stability properties 
among their peers of the same order. 
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2p 
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1 - 
0.5 
0- 
-0.5 
-1 - 
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(a) Stability region S of the im- (b) Stability region S of the ex- (c) Constrained stability region 
plicit method 



plicit method 



5„ (47) for a = 90° 



Figure 3: Stability regions for the IMEX-DIMSIM-3B pair of schemes 
6.1 Van der Pol equation 

We consider the nonhnear van der Pol equation with a split right hand side 



y 

z' 



= f{y,z)+g{y,z 
on the time interval [0, 0.5], with initial values 
y{0) = 2, z{0) = 



z 







+ 







{{l-y^)z-y)/e 



(53) 



292 



2 1,0 _ 

3 ^ ~81^ ^ 2187 



e 



1814 
19683 



(54) 



We consider e = 10 ®, a stiff case in which many methods suffer from order reduction [20 



The initialization (38) was done using the analytic derivatives. The reference solution is 
obtained with Radau-5, a stiffly accurate method [T7j, with very tight tolerances of atol = 
rtol = 5 X 10"^^. We compare the new methods with IMEX DIRK(3,4, 3), a L-stable, 
three-stage, third-order IMEX Runge-Kutta method proposed in |i4j. 

Figure |4] shows the global error, measured in the L2 norm, against step size h. A geometric 
sequence of step sizes, r, r/2, r/4 and so on, were used. Order reduction can be clearly 
observed for the IMEX Runge-Kutta method, which yields second-order convergence. The 
IMEX DIMSIM converges at the theoretical third order and gives more accurate result than 
the IMEX Runge-Kutta method. Second-order IMEX DIMSIMs also produced no order 
reduction; detailed results have been reported in |1^. These results indicate that the high 
stage order of IMEX DIMSIMs make them particularly attractive for solving stiff problems, 
where Runge-Kutta methods may suffer from order reduction. 
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Figure 4: Convergence results for third-order IMEX schemes on the van der Pol equation. 
6.2 Gravity waves 

Consider the dynamics of gravity waves, which is governed by the compressible Euler equa- 
tion in the conservative form [21] 



dp 
di 



+ V ■ (pu) 







dpu 



+ V ■ (puu + pi) 



(55a) 



dpd 
~dt 



V ■ (p^u) 



0, 



where p is the density, (u is the velocity, 9 is the potential temperature, and I is a 2 x 2 
identity matrix. The prognostic variables are (p, (pu)-^, (pO)'^)'^. The pressure p in the 
momentum equation is computed by the equation of state 



P = Po 



p9R 
Po 



(55b) 



To maintain the hydrostatic state, we follow the splitting introduced in [21] 

p(x,t) = p(z)+p'(x,t) 
(p^)(x,t) = p)(^) + (p^^)'(x,t) 
p(x,t) = p{z)+p'{x,t) 

where the reference (overlined) values are in hydrostatic balance. The gravity wave equation 



(55) can be rewritten as 



~dt 
dpu 

~dt 

djpoy 

dt 



-V ■ (pu) 

-V ■ (puu + p'l) - p'ge^ 
-V ■ (p^u) , 
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(56a) 



closed by the equation of state 



p' = Po[j^y^ -p. (56b) 

The 2D mesh is generated by the software GMSH [22]. The spatial discretization uses 
discontinuous Galerkin finite elements and was developed by Blaise et al. [19]. Figure [s] 
shows the density, velocity, potential temperature, and pressure variables after 900 seconds 
of simulation time. 

The advantage of implicit-explicit time-stepping over explicit time-stepping schemes for 
this problem has been demonstrated in [23] . To apply IMEX integration the right-hand side 



of (|56a|) is additively split into linear and nonlinear parts. The linear term 

(57) 



V ■ (pu) 

■ (p'l) + p'g^z 
V • {pen) 



with the pressure linearized as 

P = ^ (P^) 

pe 

is solved implicitly, while the remaining (nonlinear) terms are solved explicitly. 

All the experiments are performed on a workstation with 4 Intel Xeon E5-2630 Processors 
(24 cores in total) using 12 MPI threads. Note that the parallelization is not implemented at 
time-stepping level but at the spatial discretization level, therefore the parallel performance 
does not be affect the comparison of various time integrators. 

Here we compare the performance of IMEX methods for a simulation window of 30 
seconds. The second order methods are IMEX-DIMSIM-2B and L-stable, two-stage, second- 
order IMEX DIRK(2,3,2) [1]. The third order methods are IMEX-DIMSIM-3B and IMEX 
DIRK(3,4, 3) |t4j. The integrated L2 errors for all prognostic variables are measured against 
a reference solution. The reference solution was obtained by applying an explicit RK method 
to solve the original (non-split) model with a very small time step h = 0.005. 

The error versus computational effort diagrams are shown in Figure [6j All the methods 
display the theoretical orders of convergence. IMEX DIMSIMs and IMEX RK methods 
perform similarly, with IMEX DIMSIMs yielding slightly better accuracy when the same 
time steps are chosen. Also, IMEX DIMSIMs are slightly more efficient in terms of CPU time 
than the IMEX RK methods of the same order. Note that this specific DG implementation 
requires the solution to be recovered at each time step, therefore the termination procedure 
has been applied after each each time step. The implementation can be optimized such as 
to apply the termination procedure only once at the end of the simulation; this would result 
in additional savings in computational cost. As the order increases, the number of stages 
required by an IMEX RK method grows rapidly due to order conditions, while an IMEX 
DIMSIM typically uses a number of stages equal to its order. Consequently, we expect that 
IMEX DIMSIM methods will become even more competitive for higher orders. 
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uv Magnitude 




-0.0Q1538 -0.16318 

(c) potential temperature (d) pressure 

Figure 5: Solution of the gravity waves after 900 simulation seconds. The results are ob- 
tained with a third-order discontinuous Galerking space discretization and third-order IMEX 
DIMSIM time integration. 




Figure 6: Integrated L2 errors against time steps (a) and CPU time (b) for difference IMEX 
schemes. The errors are computed after 30 s of simulation. A geometric sequence of step 
sizes, r, r/2, r/4 and so on, is used. 
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7 Conclusions and future work 



In this paper introduce a new family of partitioned time integration methods based on high 
stage order general linear methods. We prove that the general linear framework is well 
suited for the construction of multi-methods. Specifically, owing to the high stage orders, no 
coupling conditions are needed to ensure the order of accuracy of the partitioned GLM. 

We apply the partitioned general linear framework to construct new implicit-explicit 
GLM pairs, together with appropriate starting and ending procedures. The linear stability 
analysis proposes the use of constrained stability functions to quantify the joint stability 
of the IMEX pair. A Prothcro-Robinson convergence analysis reveals that the order of an 
IMEX GLM scheme on very stiff problems is dictated by the stage order of its non-stiff 
component; in particular, no order reduction appears if the explicit method has a full stage 
order. This result indicates that IMEX GLMs are particularly attractive for solving stiff 
problems, where other multistage methods may suffer from order reduction. 

We discuss the construction of practical IMEX GLM pairs starting from known implicit 
schemes and adding an appropriate explicit counterpart. This strategy is applied to build 
second and third order IMEX diagonally-implicit-explicit multi-stage integration methods. 
Numerical experiments with the van der Pol equation confirm the fact that IMEX GLMs 
converge at full order while IMEX RK methods suffer from order reduction. The two dimen- 
sional gravity wave system is an important step towards solving real PDE-based problems. 
The new IMEX-DIMSIM schemes perform slightly better than the IMEX RK methods of 
the same order. 

Future work will develop IMEX-GLMs of higher orders, will endow them with adaptive 
time stepping capabilities, and will study their advantages compared to other existing IMEX 
familiess. There are also implementation issues that deserve further exploration. 
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